en 

o 

(N 



S 






Noname manuscript No. 

(will be inserted by the editor) 



Solving Heat Conduction Problems by the Direct 
Meshless Local Petrov-Galerkin (DMLPG) method 

Davoud Mirzaei • Robert Schaback 



Received: date / Accepted; date 

Abstract As an improvement of the Meshless Local Petrov-Galerkin (MLPG), 
the Direct Meshless Local Petrov-Galerkin (DMLPG) method is applied here 



C^ to the numerical solution of transient heat conduction problem. The new tech- 

'^'n nique is based on direct recoveries of test functionals (local weak forms) from 

\G values at nodes without any detour via classical moving least squares (MLS) 

^ shape functions. This leads to an absolutely cheaper scheme where the numer- 

■^r ical integrations will be done over low-degree polynomials rather than com- 

1/" . plicated MLS shape functions. This eliminates the main disadvantage of MLS 

r1 based methods in comparison with finite element methods (FEM) , namely the 

^ costs of numerical integration. 



Keywords Generalized Moving least squares (GMLS) approximation • 
Meshless methods • MLPG methods • DMLPG methods • Heat conduction 
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OO 1 Introduction 

'^ Meshless methods have received much attention in recent decades as new tools 

O 

r/-\ to overcome the difficulties of mesh generation and mesh refinement in classical 
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mesh-based methods such as the finite element method (FEM) and the finite 
volume method (FVM). 

The classification of numerical methods for solving PDEs should always 
start from the classification of PDE problems themselves into strong, weak, 
or local weak forms. The first is the standard pointwise formulation of differ- 
ential equations and boundary conditions, the second is the usual weak form 
dominating all FEM techniques, while the third form splits the integrals of 
the usual global weak form into local integrals over many small subdomains, 
performing the integration by parts on each local integral. Local weak forms 
are the basis of all variations of the Meshless Local Petrov-Galerkin technique 
(MLPG) of S.N. Atluri and collaborators 1_. This classification is dependent 
on the PDE problem itself, and independent of numerical methods and the 
trial spaces used. Note that these three formulations of the "same" PDE and 
boundary conditions lead to three essentially different mathematical problems 
that cannot be identified and need a different mathematical analysis with re- 
spect to existence, uniqueness, and stability of solutions. 

Meshless trial spaces mainly come via Moving Least Squares or kernels like 
Radial Basis Functions. They can consist of global or local functions, but they 
should always parametrize their trial functions ^^entirely in terms of nodes^^ [3 
[T5] and require no triangulation or meshing. 

A third classification of PDE methods addresses where the discretization 
lives. Domain type techniques work in the full global domain, while boundary 
type methods work with exact solutions of the PDE and just have to care for 
boundary conditions. This is independent of the other two classifications. 

Consequently, the literature should confine the term "meshless" to be a 
feature of trial spaces, not of PDE problems and their various formulations. 
But many authors reserve the term truly meshless for meshless methods that 
either do not require any discretization with a background mesh for calculating 
integrals or do not require integration at all. These techniques have a great 
advantage in computational efficiency, because numerical integration is the 
most time-consuming part in all numerical methods based on local or global 
weak forms. This paper focuses on a truly meshless method in this sense. 

Most of the methods for solving PDEs in global weak form, such as the 
Element-Free Galerkin (EFG) method [1], are not truly meshless because a 
triangulation is still required for numerical integration. The Meshless Local 
Petrov-Galerkin (MLPG) method solves PDEs in local weak form and uses no 
global background mesh to evaluate integrals because everything breaks down 
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to some regular, well-shaped and independent sub-domains. Thus the MLPG 
is known as a truly meshless method. 

We now focus on meshless methods using Moving Least Squares as trial 
functions. If they solve PDEs in global or local weak form, they still suffer from 
the cost of numerical integration. In these methods, numerical integrations are 
traditionally done over MLS shape functions and their derivatives. Such shape 
functions are complicated and have no closed form. To get accurate results, 
numerical quadratures with many integration points are required. Thus the 
MLS subroutines must be called very often, leading to high computational 
costs. In contrast to this, the stiffness matrix in finite element methods (FEMs) 
is constructed by integrating over polynomial basis functions which are much 
cheaper to evaluate. This relaxes the cost of numerical integrations. For an 
account of the importance of numerical integration within meshless methods, 
we refer the reader to [2] . 

To overcome this shortage within the MLPG based on MLS, Mirzaei and 
Schaback [8i proposed a new technique. Direct Meshless Local Petrov-Galerkin 
(DMLPG) method, which avoids integration over MLS shape functions in 
MLPG and replaces it by the much cheaper integration over polynomials. It 
ignores shape functions completely. Altogether, the method is simpler, faster 
and often more accurate than the original MLPG method. DMLPG uses a gen- 
eralized MLS (GMLS) method of |^ which directly approximates boundary 
conditions and local weak forms as some Junctionals, shifting the numerical 
integration into the MLS itself, rather than into an outside loop over calls to 
MLS routines. Thus the concept of GMLS must be outlined first in Section [2] 
before we can go over to the DMLPG in Section |4] and numerical results for 
heat conduction problems in Section [6J 

The analysis of heat conduction problems is important in engineering and 
applied mathematics. Analytical solutions of heat equations are restricted to 
some special cases, simple geometries and specific boundary conditions. Hence, 
numerical methods are unavoidable. Finite element methods, finite volume 
methods, and finite difference methods have been well applied to transient 
heat analysis over the past few decades 5]. MLPG methods were also devel- 
oped for heat transfer problems in many cases. For instance, J. Sladek et.al. 
[H] proposed MLPG4 for transient heat conduction analysis in functionally 
graded materials (FGMs) using Laplace transform techniques. V. Sladek et.al. 
[T5] developed a local boundary integral method for transient heat conduction 
in anisotropic and functionally graded media. Both authors and their collab- 
orators employed MLPG5 to analyze the heat conduction in FGMs [TOl[TT|. 
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The aim of this paper is the development of DMLPG methods for heat 
conduction problems. This is the first time where DMLPG is applied to a time- 
dependent problem. Moreover, compared to [5], we will discuss all DMLPG 
methods, go into more details and provide explicit formulae for the numerical 
implementation. DMLPGl/2/4/5 will be proposed, and the reason of ignoring 
DMLPG3/6 will be discussed. The new methods will be compared with the 
original MLPG methods in a test problem, and then a problem in FGMs will 
be treated by DMLPGl. 

In all application cases, the DMLPG method turned out to be superior to 
the standard MLPG technique, and it provides excellent accuracy at low cost. 



2 Meshless methods and GMLS approximation 

Whatever the given PDE problem is and how it is discretized, we have to find 
a function u such that M linear equations 

\k{u) = /3fc, 1 s^ fc s^ M, (1) 

defined by M linear Junctionals Xi,. . . ,Xm and M prescribed real values 
/3i, . . . , Pm are to be satisfied. Note that weak formulations will involve func- 
tionals that integrate u or a derivative against some test function. The func- 
tional can discretize either the differential equation or some boundary condi- 
tion. 

Now meshless methods construct solutions from a trial space whose func- 
tions are parametrized ^^ entirely in terms of nodes" [3]. We let these nodes form 
a set X := {xi, . . . ,X]\[}. Theoretically, meshless trial functions can then be 
written as linear combinations of shape functions ui , . . . , un with or without 
the Lagrange conditions Uj{xk) — Sjk, 1 ^ j,k < N a.s 

N 



uix) = ^Uj(a;)u(a;j) 



in terms of values at nodes, and this leads to solving the system (fTl) in the 
form 

N 

^k{u) ^^\k{uj)u{xj) = /3fc, l^ki^M 

approximately for the nodal values. Setting up the coefficient matrix requires 
the evaluation of all functionals on all shape functions, and this is a tedious 
procedure if the shape functions are not cheap to evaluate, and it is even 
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more tedious if the functionals consist of integrations of derivatives against 
test functions. 

But it is by no means mandatory to use shape functions at this stage at 
aU. If each functional A^ can be well approximated by a formula 

N 

Afc(u) « ^ ajku{xj) (2) 

in terms of nodal values for smooth functions u, the system to be solved is 

]ajkuix^)^Pk, l^k^M (3) 



N 



i=i 



without any use of shape functions. There is no trial space, but everything is 
still written in terms of values at nodes. Once the approximate values u{xj) at 
nodes are obtained, any multivariate interpolation or approximation method 
can be used to generate approximate values at other locations. This is a post- 
processing step, independent of PDE solving. 

This calls for efficient ways to handle the approximations (pi) to functionals 
in terms of nodal values. We employ a generalized version of Moving Least 
Squares (MLS), adapted from [5|, and without using shape functions. 

The techniques of [51 and [51 allow to calculate cocfhcients ajk for ^ very 
effectively as follows. We fix k and consider just A := Afe. Furthermore, the set 
X will be formally replaced by a much smaller subset that consists only of the 
nodes that are locally necessary to calculate a good approximation of A^, but 
we shall keep X and N in the notation. This reduction of the node set for the 
approximation of Xk will ensure sparsity of the final coefficient matrix in pi). 

Now we have to calculate a coefficient vector a{Xk) = {cnik, ■ ■ ■ , ctNk)'^ € 
M^ for ([2]) in case of A = Afc. We choose a space V of polynomials which is 
large enough to let zero be the only polynomial p in V that vanishes on X. 
Consequently, the dimension Q oi V satisfies Q ^ N, and the Q x N matrix 
P of values Pi{xj) of a basis pi, ■ ■ ■ ,pq of V has rank Q. Then for any vector 
w = {wi, . . . , wn)^ of positive weights, the generalized MLS solution a(A) to 
(pi) can be written as 

a(Xk) - WP^{PWP^)-^Xk{V) (4) 

where W is the diagonal matrix with diagonal w and Xk{V) G MP is the vector 
with values Xk{pi) , ■ ■ ■ , Xkipq) ■ 

Thus it suffices to evaluate A^ on low-order polynomials, and since the 
coefficient matrix in ([4| is independent of fc, one can use the same matrix 
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for different Afe as long as X does not change locally. This will significantly 
speed up numerical calculations, if the functional A^ is complicated, e.g. a 
numerical integration against a test function. Note that the MLS is just behind 
the scene, no shape functions occur. But the weights will be defined locally 
in the same way as in the usual MLS, e.g. we choose a continuous function 
<j) : [0, cx)) -> [0, oo) with 

- 0(r) >0, 0<r < 1, 

- 0(r) = 0, r ^ 1, 

and define 

Wj^x) — 4> 



Cj||2 



5 
for (5 > as a weight function, if we work locally near a point x. 

3 MLPG Formulation of Heat Conduction 

In the Cartesian coordinate system, the transient temperature field in a het- 
erogeneous isotropic medium is governed by the diffusion equation 

p{x)c{x) — (x, i) = V • {nVu) + fix, t), (5) 

where x € f2 and ^ t ^ tp denote the space and time variables, respectively, 
and tp is the final time. The initial and boundary conditions are 



w(x,0) = 


= uo{x), X £ n, 




(6) 


u{x,t) = 


= UD{x,t), XETd, 


Os^ti^tF, 


(7) 


Kix)-^ix,t) -- 


^UN{x,t), xeTN, 


< t 5$ tF- 


(8) 



In (l5])-(|8|, u{x,t) is the temperature field, k{x) is the thermal conductivity 
dependent on the spatial variable x, p{x) is the mass density and c{x) is the 
specific heat, and /(x, t) stands for the internal heat source generated per unit 
volume. Moreover, n is the unit outward normal to the boundary F , ujj and 
UN are specified values on the Dirichlet boundary Fd and Neumann boundary 
Iat where F = F^ U F^- 

Meshless methods write everything entirely in terms of scattered nodes 
forming a set X = {xi,X2, ■ ■ ■ ,xis[} located in the spatial domain J7 and 
its boundary F. In the standard MLPG, around each Xk a small subdomain 
J7^ C J7 = i7 U F is chosen such that integrations over Q^ are comparatively 
cheap. For instance, Q^ is conveniently taken to be the intersection of Q with 
a ball B{xk, vq) of radius rg or a cube (or a square in 2D) S{xk, tq) centered 
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at Xk with side- length tq. On these subdoniains, the PDE including boundary 
conditions is stated in a localized weak form 



d_ 

dt 



/ pcuvdf2= / V-(KVu)udi7+ / fvdf2, (9) 



for an appropriate test function v. Applying integration by parts, this weak 
equation can be partially symmetrized to become the first local weak form 

— / pcuvdQ^ / K-^vdr- / nVu-Vvdn + / fvdn. (10) 

The second local weak form, after rearrangement of ([5]) and integration by 
parts twice, can be obtained as 

d f 1 f f dv f du 

-^r- / —pcuvdfi— / uAv dfi — / u-;— dr + I v-^^dF 
ot Jnk K jQk jQQk On jQQk dn 

/■ 1 /■ 1 ^^^' 

+ / -VK-Vuvdl2+ / -fvdQ. 
Jnk K J Qk K 

If the boundary of the local domain il^ hits the boundary of /?, the MLPG 
inserts boundary data at the appropriate places in order to care for boundary 
conditions. Since these local weak equations are all affine-linear in u even 
after insertion of boundary data, the equations of MLPG are all of the form 
(fTj) after some rearrangement, employing certain linear functionals A^;. In all 
cases, the MLPG evaluates these functionals on shape functions, while our 
DMLPG method will use the GMLS approximation of Section [5] without any 
shape function. 

However, different choices of test functions v lead to the six different well- 
known types of MLPG. The variants MLPGl/5/6 are based on the weak 



formulation (10 1. If u is chosen such that the first integral in the right hand 
side of (10) vanishes, we have MLPGl. In this case v should vanish on df2^. 
If the Heaviside step function v on local domains is used as test function, the 
second integral disappears and we have a pure local boundary integral form in 
the right hand side. This is MLPG5. In MLPG6, the trial and test functions 
come from the same space. MLPG2/3 are based on the local unsymmetric weak 
formulation (|9|. MLPG2 employs Dirac's delta function as the test function in 
each n^, which leads to a pure collocation method. MLPG3 employs the error 
function as the test function in each J7^. In this method, the test functions can 
be the same as for the discrete least squares method. The test functions and 
the trial functions come from the same space in MLPG3. Finally, MLPG4 (or 



LBIE) is based on the weak form (111, and a modified fundamental solution 
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of the corresponding elliptic spatial equation is employed as a test function in 
each subdomain. 

We describe these types in more detail later, along with the way we modify 
them when going from MLPG to DMLPG. 



4 DMLPG Formulations 

Independent of which variation of MLPG we go for, the DMLPG has its special 
ways to handle boundary conditions, and we describe these first. 

Neither Lagrange multipliers nor penalty parameters are introduced into 
the local weak forms, because the Dirichlet boundary conditions are imposed 
directly. For nodes Xk G r^, the values u(xk,t) ~ U£i{xk,t) are known from 
the Dirichlet boundary conditions. To connect them properly to nodal values 
u{xj,t) in neighboring points Xj inside the domain or on the Neumann bound- 
ary, we turn the GMLS philosophy upside down and ask for coefficients aj{xk) 
that allow to reconstruct nodal values at x^. from nodal values at the Xj. This 
amounts to setting X^ — S^,, in Section [2] , and we get localized equations for 
Dirichlet boundary points Xk as 

N 

^aj{xk)u{xj,t) =UDixk,t), Xk^Fn, ie[0,iF]- (12) 



Note that the coefficients are time-independent. In matrix form, (12 1 can be 
written as 

Bu{t)^UD{t), (13) 

where u{t) £ M.-^ is the time-dependent vector of nodal values at xi, X2, ■■■, xn- 
These equations are added into the full matrix setup at the appropriate places, 
and they are in truly meshless form, since they involve only values at nodes 



and are without numerical integration. Note that ( 10 1 has no integrals over 



the Dirichlet boundary, and thus we can impose Dirichlet conditions always 



in the above strong form. For (11) there are two possibilities. We can impose 



the Dirichlet boundary conditions either in the local weak form or in the 



collocation form (12 1. Of course the latter is the cheaper one. 

We now turn to Neumann boundary conditions. They can be imposed in 
the same way as Dirichlet boundary conditions by assuming Afe(u) ~ ^i^k) 
in the GMLS approximation 

'^a-i{xk)u{xj,t) = ---{xk,t), XkETN, te[0,i_F]. (14) 



DMLPG for heat conduction 



Note that the coefRcients again are time-independent, and we get a Unear 



system Hke ( 13 1, but with a vector Uj^{t) of nodal values of normal derivatives 
in the right-hand side. This is collocation as in subsection |4.2[ But it is often 
more accurate to impose Neumann conditions directly into the local weak 



forms (10 1 and (11). We will describe this in more detail in the following 
subsections. We now turn the different variations of the MLPG method into 
variations of the DLMPG. 



4.1 DMLPGl/5 



These methods are based on the local weak form (10). This form recasts to 



— / pcuv dil + / KWu-Wvdil— / K—-vdr 

n (15) 

u^v dr + fv dfl 
rjvn9r2j Jfij 

after inserting the Neumann boundary data from (p|, when the domain Q^ of 



(10) hits the Neumann boundary Fm- All integrals in the top part of ( |15[ ) can 
be efficiently approximated by GMLS approximation of Section [2] as purely 
spatial formulas 



Ai.fclw) := / pcuvdil Ri Xi,k{u) = y^ ai^j{xk)u{xj), 

^2,k{u) -.^ ~ KVu-VvdQ ~ X2.k{u) ^^S2o,2.]{xk)u{xj), (16) 

•^''' 3 = 1 

r d — -— - ^ 

A3.fc(w) ■■= ~ K—-vdr sa A3.fe(u) = ^ a:i,j{xk)u{xj). 

Jdni\r^ On ' ~[ ' 

While the two others can always be summed up, the first formula, if applied 
to time-varying functions, has to be modified into 

d " 



r d 

/ pcuvdf2 KiS^ aij{xk)Tr:u{xj,t) 
Jni' ~f ot 



and expresses the main PDE term not in terms of values at nodes, but rather 
in terms of time derivatives of values at nodes. 

Again, everything is expressed in terms of values at nodes, and the coeffi- 
cients are time-independent. Furthermore, Section [2] shows that the u part of 
the integration runs over low-order polynomials, not over any shape functions. 
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The third functional can be omitted if the test function v vanishes on 
dil^ \ IV- This is DMLPGl. An example of such a test function is 

'\\x~ Xkh^ 



V = v{x;xk) = 



ro 



where cf) is the weight function in the MLS approximation with the radius 6 of 
the support of the weight function being replaced by the radius ro of the local 
domain il^. 

In DMLPG5, the local test function is the constant v = 1. Thus the func- 



tionals A2,fe of (16) are not needed, and the integrals for Ai^fc take a simple 
form, if c and p are simple. DMLPG5 is slightly cheaper than DMLPGl, be- 
cause the domain integrals of A2,fe are replaced by the boundary integrals of 

As^fc- 

Depending on which parts of the functionals are present or not, we finally 
get a time-dependent system of the form 



A(i)|-M(i)+v4(^)M(i) = (i), ^ = 2 or 3 



(17) 



where u(t) is the time-dependent vector 

u{t) = (u(xi, i), . . . , u{xN, t)f e M^ 

of nodal values, {t) £ R^^ collects the time-dependent right-hand sides with 
components 

bk 

and A[^^ = ag^jixk), i^l,2, 3. The fc-th row of A(^) is 

ai'^ ^WPiPWP'^y'Xi^kiV), £=1,2,3, 



f{x,t)v{x]Xk)dn + i UN{x,t)v{x;xk)dr, 



pcpiv dil, / pcp2V dfi, . . . , / 
r2j J n'; Jfi 



H T 



ni 



PCPQ' 



!d[2 



where 

A2,fc(P) 

As we can immediately see, numerical integrations are done over low-degree 
polynomials pi,p2, ■■■,Pq only, and no shape function is needed at all. This 
reduces the cost of numerical integration in MLPG methods significantly. 



/ nVpi ■ Vw dil, / k\'p2 ■ Vw df2, . . . , / nVpq ■ Vw dil 
J f}^ J f2^ J u^ 

,^,dP ( -^^^^-^ ( -^^^ 



dn 



K — — V dP, 
on 



-\ T 



K^-^vdP 
on 
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4.2 DMLPG2 

In this method, the test function v on the local domain Q^ in (|9| is replaced 
by the test functional Sx^ , i.e. we have strong collocation of the PDE and all 
boundary conditions. Depending on where Xk lies, one can have the functionals 

M2,fe(w) := ^-(pcM)(xfe), (18) 

on 

m,k{u) := V • {KVu){xk) 
connecting u to Dirichlet, Neumann, or PDE data. The first form is used on 



the Dirichlet boundary, and leads to (12 1 and (13). The second applies to 



points on the Neumann boundary and is handled by ( 14 ) , while the third can 
occur anywhere in Q independent of the other possibilities. In all cases, the 
GMLS method of Section [2] leads to approximations of the form 

N 

M»,fe(") ~ '^a,i^j{xk)u{x^j), i = 1,2, 3 

entirely in terms of nodes, where values on nodes on the Dirichlet boundary 
can be replaced by given data. 

This DMLPG2 technique is a pure collocation method and requires no 
numerical integration at all. Hence it is truly meshless and the cheapest among 
all versions of DMLPG and MLPG. But it needs higher order derivatives, and 
thus the order of convergence is reduced by the order of the derivative taken. 
Sometimes DMLPG2 is called Direct MLS Collocation (DMLSC) method [Hj. 

It is worthy to note that the recovery of a functional such as H2,k{u) or 



fJ-a.kiu) in (18) using GMLS approximation gives GMLS derivative approxima- 
tion. These kinds of derivatives have been comprehensively investigated in |5] 
and a rigorous error bound was derived for them. Sometimes they are called 
diffuse or uncertain derivatives, because they are not derivatives of shape func- 
tions, but [9] proves there is nothing diffuse or uncertain about them and they 
are direct and usually very good numerical approximation of corresponding 
function derivatives. 



4.3 DMLPG4 



This method is based on the local weak form (11) and uses the fundamental 
solution of the corresponding elliptic spatial equation as test function. Here we 
describe it for a two-dimensional problem. To reduce the unknown quantities 
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in local weak forms, the concept of companion solutions was introduced in |16j . 
The companion solution of a 2D Laplace operator is 

v{x;y) = --\n—, r ^ \\x - y\\2, 
ZTT r 

which corresponds to the Poisson equation Av{x; y) + 5{r) = and thus is 
a fundamental solution vanishing for r — r^. Dirichlet boundary conditions 



for DMLPG4 are imposed as in (12 1. The resulting local integral equation 



corresponding to a node Xk located inside the domain or on the Neumann 
part of the boundary is 

d f 1 ,„ , , f dv 



1 f dv f I 

. —pcuv dQ—aku{xk) + -f- ^—udF— / —Vn-Vuvdfl 
at J Ok K JdO'' on J„k K 



f 1 

UNvdr+ / —Jvdfl, 



(19) 



where ak is a coefficient that depends on where the source point Xk lies. It is 
1/2 on the smooth boundary, and 0fe/(27r) at a corner where the interior angle 
at the point Xk is 0k- The symbol j^ represents the Cauchy principal value 
(CPV). For interior points Xk we have ak — 1 and CPV integrals are replaced 
by regular integrals. 
In this case 



Ai,fe(^) 



1 /■ 1 

—pcpivdQ, / —pcp2vdf2, 



1 



-\T 



pcpQV df2 



,W 



and X2MV) = o^kXmT) + X'^HV) + X^V), where 



(2), 



.(3)/ 



^2,kiT^) = [PliXk),P2{Xk), ■ ■ ■ ,PQiXk)y 






dv 



Pi dr, f --P2 dr,...,f 

rk On Jpk on Jf 



n 



dv 
dn 



PQ 



dr 



1 /■ 1 

—S/K-Vpivdn, / —'s/K-Vp2vdn, 

fife ^ Jn'' '* 



m 1^ 



-Vk • VpQ V dn 



Finally, we have the time-dependent linear system of equations 

A(l)|«(i)+A(2Vi) = (t), 

where the fc-th row of A^^^ is 



(20) 



An 



wPiPWP'y'XeA'P), ^ = 1,2. 
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The components of the right-hand side are 

bk{t) = ——-f{x,t)v(x]Xk)dQ+ I UN{x,t)v{x;xk)dr. 

jQk K(x) JafijnFjv 

This technique leads to weakly singular integrals which must be evaluated by 
special numerical quadratures. 



4.4 DMLPG3/6 

In both MLPG3 and MLPG6, the trial and test functions come from the 
same space. Therefore they are Galerkin type techniques and should better 
be called MLG3 and MLG6. But they annihilate the advantages of DMLPG 
methods with respect to numerical integration, because the integrands include 
shape functions. Thus we ignore DMLPG3/6 in favour of keeping all benefits 
of DMLPG methods. Note that MLPG3/6 are also rarely used in comparison 
to the other MLPG methods. 



5 Time Stepping 

To deal with the time variable in meshless methods, some standard meth- 
ods were proposed in the literature. The Laplace transform method [T^llOj . 
conventional finite difference methods such as forward, central and backward 
difference schemes are such techniques. A method which employs the MLS 
approximation in both time and space domains, is another different scheme 

In our case the linear system (l3| turns into the time-dependent version 



(17) coupled with (13) that could, for instance, be solved like any other linear 
first-order implicit Differential Algebraic Equations (DAE) system. Invoking 
an ODE solver on it would be an instance of the Method of Lines. If a conven- 
tional time-difference scheme such as a Grank-Nicolson method is employed, if 
the time step At remains unchanged, and if M = N, then a single LU decom- 
position of the final stiffness matrix and corresponding backward and forward 
substitutions can be calculated once and for all, and then the final solution 
vector at the nodes is obtained by a simple matrix-vector iteration. 

The classical MLS approximation can be used as a postprocessing step to 
obtain the solution at any other point x d f2. 
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6 Numerical results 

Implementation is done using the basis polynomials 



{x - zf 
h\0\ 



0^|/3Km 



where h is an average mesh-size, and z is a fixed evaluation point such as a 
test point or a Gaussian point for integration in weak-form based techniques. 
Here /3 = (/3i, . . . , ^d) G N^ is a multi-index and |/3| = /3i + . . . + Pa- If 
X — (xi, • • ■ , Xd) then x'^ — Xi ■ ■ ■ Xd ■ T^^^s choice of basis function, instead 
of {a;^}o^|;3|^mi leads to a well-conditioned matrix PWP^ in the (G)MLS 
approximation. The effect of this variation on the conditioning has been ana- 
lytically investigated in j9]. 

A test problem is first considered to compare the results of MLPG and 
DMLPG methods. Then a heat conduction problem in functionally graded 
materials (FGM) for a finite strip with an exponential spatial variation of ma- 
terial parameters is investigated. In numerical results, we use the quadratic 
shifted scaled basis polynomial functions (m = 2) in (G)MLS approximation 
for both MLPG and DMLPG methods. Moreover, the Gaussian weight func- 
tion 



||x- a;j||2 > 5 



where 6 ~ S^h and c = cg/i is used. The parameter 5q should be large enough to 
ensure the regularity of the moment matrix PWP"^ in (G)MLS approximation. 
It depends on the degree of polynomials in use. Here we put 5q = 2m. The 
constant cg controls the shape of the weight function and has influence on the 
stability and accuracy of (G)MLS approximation. There is no optimal value 
for this parameter at hand. Experiments show that 0.4 < cq < 1 lead to more 
accurate results. 

All routines were written using Matlab® and run on a Pentium 4 PC 
with 2.50 GB of Memory and a twin-core 2.00 GHz CPU. 
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6.1 Test problem 

Let n = [0, 1]^ C M^ and consider Equations ([5])-([8| with pc = 27r^, k = 1 and 
f{x,t) = 0. Boundary conditions using x — (X11X2) ^ ^"^ ^re 

du , , , , , , 

^ = 0, (xi, X2 = 0) U (xi, X2 = 1), XI G [0, 1], 

u = e"*cos(7rx2), (xi =0,X2), X2 e [0,1], 
u = -e"* cos(7rx2), (xi = 1, X2), X2 G [0, 1]- 

The initial condition is u(x, 0) — cos(7rxi) cos(7rx2), and u{x, t) — e^* cos(7rxi) cos(7rx2) 
is the exact solution. Let tp = 1 and At — 0.01 in the Crank-Nicolson scheme. 
A regular node distribution with distance h in both directions is used. In 
Table[l]the CPU times used by MLPGl/2/4/5 and DMLPGl/2/4/5 are com- 
pared. As we can immediately see, DMLPG methods are absolutely faster 
than MLPG methods. There is no significant difference between MLPG2 and 
DMLPG2, because they are both collocation techniques and no numerical in- 
tegration is required. 

INSERT TABLE 1 

The maximum absolute errors are drawn in FigurefTland compared. MLPG2 
and DMLPG2 coincide, but DMLPGl/4/5 are more accurate than MLPGl/4/5. 
DMLPGl is the most accurate method among all. Justification needs a rig- 
orous error and stability analysis which is not presented here. But, accord- 
ing to [SlIS] and all numerical results, we can expect an error behavior like 
^j-^m+i-fe-j^ where k is the maximal order of derivatives of u involved in 
the functional, and if numerical integration and time discretization have even 
smaller errors. 

INSERT FIGURE 1 

For more details see the elliptic problems in [8; where the ratios of errors of 
both method types are compared for m — 2,3, 4. 



6.2 A problem in FGMs 

Consider a finite strip with a unidirectional variation of the thermal conduc- 
tivity. The exponential spatial variation is taken 

K(a;) = Koexp(7xi), (21) 
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with Ko = 17 Wm^^ °C^^ and pc ~ 10^. This problem has been considered 
in [T^ using the meshlcss LBIE method (MLPG4) with Laplace transform in 
time, and in [71IS] using MLPG4/5 with MLS approximation for both time and 
space domains, and in [T3] using a RBF based meshless collocation method 
with time difference approximation. 

In numerical calculations, a square with a side-length a ~ A cm and a 
11 X 11 regular node distribution is used. 

Boundary conditions are imposed as bellow: the left side is kept to zero 
temperature and the right side has the Heaviside step time variation i.e., u — 
TH{t) with T = 1°C . On the top and bottom sides the heat flux vanishes. 

We employed the ODE solver ode 15s from MATLAB for the final DAE 
system, and we used the relative and absolute tolerances le-5 and le-6, re- 
spectively. With these, we solved on a time interval of [0 60] with initial 
condition vector mq at time 0. The Jacobian matrix can be defined in advance 
because it is constant in our linear DAE. The integrator will detect stiffness 
of the system automatically and adjust its local stepsize. 

In special case with an exponential parameter 7 = which corresponds to 
a homogeneous material the analytical solution 

Txi , 2 -^Tcosmr . niixi ( aon'^Tr'^t 

u[x, t) — ! — > sm X exp 

a TT ^-^ n a \ a^ 

n—l ^ 

is available. It can be used to check the accuracy of the present numerical 
method. 

Numerical results are computed at three locations along the xi-axis with 
Xi/a — 0.25, 0.5 and 0.75. Results are depicted in Fig. l2] An excellent agree- 
ment between numerical and analytical solutions is obtained. 

INSERT FIGURE 2 

It is known that the numerical results are rather inaccurate at very early 
time instants and at points close to the application of thermal shocks. There- 
fore in Fig. [3] we have compared the numerical and analytical solutions at very 
early time instants (t S [0, 0.4]). Besides, in Fig.Elthe numerical and analytical 
solutions at points very close to the application of thermal shocks are given 
and compared for sample time t{70) « 10.5 sec. 

INSERT FIGURES 3 AND 4 

The discussion above concerns heat conduction in homogeneous materials 
in a case where analytical solutions can be used for verification. Consider 
now the cases 7 = 0, 20, 50, and 100 m^^, respectively. The variation of 



DMLPG for heat conduction 17 

temperature with time for the three first 7- values at position xija = 0.5 are 
presented in Fig.[5J The results are in good agreement with Figure 11 presented 
in [13], Figure 6 presented in [7^ and Figure 4 presented in [5]. 

INSERT FIGURE 5 

In addition, in Fig. |6] numerical results are depicted for 7 = 100 m~^. For 
high values of 7, the steady state solution is achieved rapidly. 

INSERT FIGURE 6 

It is found from Figs. [5] and [6] that the temperature increases with an 
increase in 7- values. 

For the final steady state, an analytical solution can be obtained as 

I ^ . \ T^<^xp(-7xi) - 1 / XI ^ „ 

u(a;,t — > 00) = J -. r — , u — > J — , as 7 — >■ 

exp(— 7a) — 1 V a 

Analytical and numerical results computed at time i = 60 sec. are presented 
in Fig. [7j Numerical results are in good agreement with analytical solutions 
for the steady state temperatures. 

INSERT FIGURE 7 
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Type 1 




Type 2 




Type 4 




Type 5 




h 


MLPG 


DMLPG 


MLPG 


DMLPG 


MLPG 


DMLPG 


MLPG 


DMLPG 


0.2 


4.3 


0.2 


0.2 


0.2 


1.9 


0.2 


1.4 


0.2 


0.1 


22.6 


0.3 


0.3 


0.3 


9.8 


0.3 


6.8 


0.3 


0.05 


116.4 


1.4 


0.8 


0.6 


52.9 


1.1 


35.6 


1.2 


0.025 


855.8 


9.6 


8.3 


7.0 


446.5 


8.3 


302.2 


8.5 



Table 1 Comparison of MLPG and DMLPG methods in terms of CPU times used (Sec.) 




-5 
10 r 



--^- MLPGl 
— ►— DMLPGl 
-■Q- MLPG2 
— •— DMLPGl 



~Et- MLPG4 
-m— DMLPG4 
--^- MLPG5 
-^— DMLPG5 



Fig. 1 Comparison of MLPG and DMLPG methods in terms of maximum errors. 
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Fig. 2 Time variation of tlie temperature at three positions with 7 = 0. 
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Fig. 3 Accuracy of method for early time instants at position xi/a = 0.5. 
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Fig. 4 Accuracy of method for points close to the thermal shock at time t(70) sec. 
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Fig. 5 Time variation of the temperature at position xi/a = 0.5 for 7 = 0, 20, 50 m ^ . 
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Fig. 6 Time variation of the temperature at positions xi/a = 0.25,0.5,0.75 for 7 = 100 
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Fig. 7 Distribution of temperature along xi-axis under steady-state loading conditions. 



